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Abstract 

A  three-dimensional  numerical  simulation  of  the  transient  response  of  a  polymer  electrolyte  membrane  fuel  cell  (PEMFC)  subjected  to 
a  variable  load  is  presented.  The  model  parameters  are  typical  for  a  laboratory- scale  cell  with  a  serpentine  flow  path  and  a  10  cm2  reactive 
area.  The  simulation  uses  a  commercial  computational  fluid  dynamics  (CFD)  solver  modified  to  include  the  electrochemical  behavior.  The 
predictions  are  based  on  an  isothermal  set  of  equations  and  include  transient  responses  of  the  cell  in  terms  of  local  distributions  of  the  current 
density  and  gas  mole  fractions.  The  predictions  show  transients  in  the  current  density  that  overshoot  the  final  state  value  when  the  cell  voltage 
is  abruptly  changed  from  0.7  to  0.5  V  for  fixed  excess  initial  stoichiometric  flowrates.  The  fixed  flowrates  are  excess  because  they  correspond 
to  stoichiometries  of  2.6  and  4.4  at  0.7  V  for  the  0.35  A/cm2  predicted  initial  current  density.  The  percent  overshoot  decreases  with  the  rate  of 
voltage  change  and  it  is  shown  to  change  with  anode  gas  flow  rates.  Also  the  magnitude  of  this  overshoot  and  undershoot  can  be  adjusted  by 
changing  the  rate  of  voltage  change  and  the  operating  conditions.  The  overshoot  behavior  for  these  excess  stoichiometric  flowrates  is  shown 
to  depend  on  changes  in  the  oxygen  mole  fraction  distributions. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

The  control,  design,  and  optimum  operation  of  polymer 
electrolyte  membrane  fuel  cells  (PEMFCs)  will  require  an 
understanding  of  its  dynamics  when  there  are  changes  in 
either  current,  voltage,  or  power.  These  dynamics  would  be 
important  for  power  conditioning  in  residential  applications 
and  for  automotive  systems  performing  in  a  Federal  Urban 
Driving  Scenario.  While  some  of  these  dynamics  could  be 
determined  from  experiments,  a  model  capable  of  predicting 
the  transient  response  would  be  useful  in  the  design  flow  fields 
and  for  optimization  of  control  schemes.  We  are  particularly 
interested  in  overshoot  behavior  because  recent  experiments 
in  our  laboratory  [1-3]  indicate  the  possibility  of  overshoot 
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behavior  dependent  on  the  rate  of  voltage  change,  and  stoi¬ 
chiometry. 

As  a  guide  for  simple  one-dimensional  (ID)  and  two- 
dimensional  (2D)  models  that  may  be  developed  in  the 
future,  we  present  here  a  three-dimensional  (3D)  solution 
to  the  isothermal,  time- dependent  phenomena  that  includes 
Navier-Stokes  equations  for  the  flow  channel  and  diffusion 
layers  of  a  serpentine  flow  field  and  the  interaction  of  theses 
equations  with  anode  and  cathode  kinetics,  and  the  water  bal¬ 
ance  throughout  the  cell.  It  is  important  to  note  that  simpler 
models  (i.e.,  2D  transient)  could  not  accurately  provide  pre¬ 
dictions  of  the  serpentine  flow-field.  Commonly  use  to  test 
membrane  electrode  assemblies  (MEA).  To  obtain  this  solu¬ 
tion,  the  3D  model  of  Shimpalee  et  al.  [4,5]  was  extended 
by  including  the  time  dependent  analysis  of  a  10-cm2  reac¬ 
tive  area.  This  set  of  equations,  as  shown  below,  does  not 
account  explicitly  for  the  condensation  of  water  vapor  result¬ 
ing  from  the  local  partial  pressure  exceeding  the  saturation 
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pressure  of  water,  but  rather  it  calculates  the  activity  of  water 
in  the  gas  phase  and  allows  the  analyst  to  manually  turn  off 
the  reaction  in  the  cells  where  the  thickness  of  condensed 
liquid  restricts  mass  transfer  to  the  catalyst.  Recent  steady- 
state  versions  of  this  model  [6-8]  account  for  this  water  phase 
change  and  include  the  energy  balance  to  predict  temperature 
changes,  but  under  some  operational  conditions  liquid  water 
does  not  significantly  affect  the  performance  of  the  cell  and 
computational  time  can  be  saved  by  neglecting  these  phase 
changes.  Hence,  an  isothermal  single-phase  approximation 
is  used  here. 

The  equations  of  this  paper  and  the  results  discussed  below 
can  be  related  to  previous  works.  For  example,  Um  et  al.  [9] 
presented  similar  equations  for  their  3D  model  and  applied 
those  to  their  straight  channel  geometries.  They  used  their 
own  code  and  the  equation  solver  was  not  specified.  Recently 
Wang  and  Wang  [10]  showed  dynamic  study  on  PEMFC 
using  commercial  CFD  code.  However,  their  geometry  was 
single  straight  channel  flow  field  that  is  not  normally  used  in 
fuel  cell  application,  therefore  the  effect  of  pressure  driven 
flow  underneath  the  gas  diffusion  layers  has  been  ignored. 
In  the  results  presented  below  we  also  used  a  commercially 
available  computational  fluid  dynamics  (CFD)  solver  and  we 
had  to  supply  subroutines  to  account  for  the  electrochemi¬ 
cal  reactions  of  hydrogen  and  oxygen  for  source  terms  in  the 
transport  equations.  The  complete  equations  are  solved  with 
a  control  volume  based  discretization  of  the  computational 
domain  to  obtain  the  velocity  and  pressure  distribution  in 
the  flow  channels  and  the  gas  diffusion  layers  for  every  time 
interval. 

In  contrast  to  the  transient  analysis  presented  here  and 
recent  paper  from  reference  [10],  preceding  modeling  stud¬ 
ies  have  focused  on  the  steady  state  behavior  of  the  fuel 
cell.  Some  models  considered  only  ID  simulations  [11-13], 
several  models  concentrated  on  2D  flow  with  transport  of 
the  reactants  and  products  in  the  flow  channels  and  across 
the  membranes  [14-21],  and  recent  models  put  emphasis  on 
the  3D  simulations  [4-10].  Some  of  These  works  have  been 
reviewed  by  Shimpalee  et  al.  [22].  One  recent  ID  transient 
model  has  been  presented  and  transient  measurements  of  the 


voltage  as  the  current  was  changed  were  used  to  estimate 
parameter  of  a  ME  A.  However,  that  model  was  not  used  to 
explain  changing  in  the  overshoot  behavior  with  stoichiom¬ 
etry.  It  may  also  to  note  that  the  non-isothermal  form  of  the 
model  presented  here  has  been  verified  with  water  balance 
measurement  recently  [7]. 

2.  Model  development 

In  this  study,  the  3D  model  of  Shimpalee  and  cowork¬ 
ers  [4,5]  is  extended  by  including  the  accumulation  (time 
dependent)  term.  Thus,  this  is  a  transient,  3D,  isothermal, 
single  phase,  and  multi-species  investigation  of  a  single  PEM 
fuel  cell  with  twenty  serpentine  channels.  Again,  an  isother¬ 
mal  set  of  equations  is  used  to  increase  computational  speed 
and  to  provide  a  basis  of  comparison  for  future  work.  The 
flow  path  consists  of  a  serpentine  gas  channel  and  the  details 
of  the  computational  domain  have  been  shown  previously, 
[5,7].  Fig.  1  also  shows  the  channel  geometry  and  associ¬ 
ated  coordinate  system.  A  thin  membrane  electrode  assembly 
is  sandwiched  between  anode  and  cathode  diffusion  layers. 
Fig.  2  shows  more  details  of  the  computational  domain,  which 
consists  of  the  anode  flow  channel,  the  anode  diffusion  layer, 
the  MEA,  the  cathode  diffusion  layer,  and  the  cathode  flow 
channel.  Fig.  2  further  shows  different  z-locations  that  are 
used  in  defining  source  terms  (see  Table  1  for  a  list  of  sym¬ 
bols  and  Tables  2  and  3  for  the  equations  used  in  this  paper). 
The  species  considered  are  hydrogen,  oxygen,  nitrogen,  and 
water  vapor.  The  fuel  cell  operation  is  characterized  as  gas 
transport  and  transformation  of  one  species  to  the  other. 
The  hydrogen  from  the  anode  flow  channel  is  transported 
through  the  diffusion  layer  toward  the  membrane.  Hydro¬ 
gen  molecules  are  dissociated  to  protons  and  electrons  in  the 
catalyst.  The  water  that  impregnates  the  MEA  hydrates  the 
protons  and  it  is  transported  by  both  electro-osmosis  and  dif¬ 
fusion  according  to  Eq.  (15)  of  Table  3.  The  air  mixture  in 
the  cathode  channel  is  transported  through  the  diffusion  layer 
toward  the  membrane  where  oxygen  reacts  with  protons. 
The  water  activity  in  the  membrane  is  simulated  by  surface- 


Fig.  1.  The  picture  shows  actual  flow-field  plate  with  the  gas  channel  and  its  geometry  model.  There  are  20  straight  channels  connected  in  a  serpentine  fashion. 
Anode  and  cathode  side  flow  channels  are  symmetric  and  placed  properly  aligned  on  top  of  each  other  [2,3]. 
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Fig.  2.  The  detail  of  computational  domain  and  grid  arrangement  used  in  this  model  [2,3].  Note  the  number  of  volume  elements  used  in  references  [2,3]  and 
shown  in  this  figure  is  greater  in  the  z-direction  than  the  number  used  here  as  discussed  in  the  text. 


based  source  terms  in  the  control  volumes  in  contact  with  the 
membrane. 

Table  2  shows  the  model  equations.  The  steady  state 
aspects  of  these  equations  have  been  discussed  elsewhere 
[4,5,22].  It  may  be  important  to  note  that  the  time  depen¬ 
dent  conservation  of  mass  equation  (Table  2,  Eq.  (1))  in  the 
3D  flow  domain  is  modified  to  include  the  electrochemi¬ 
cal  reactions  of  a  fuel  cell  by  using  the  respective  source 
terms,  Sm,  specified  through  Eqs.  (7)  and  (9)— (12)  of  Table  2. 
Note  that  the  source  terms  are  zero  in  most  of  the  computa¬ 
tional  domain.  These  terms  corresponds  to  the  consumption 
of  hydrogen  in  the  anode,  and  the  consumption  of  oxygen 
and  production  of  water  in  the  cathode.  The  flux  of  water 
is  also  included  as  a  source  term  at  the  anode  and  cathode 
(i.e.,  Eqs.  (10)  and  (12))  by  accounting  for  the  diffusion, 
water  content  in  the  membrane,  and  electro- osmotic  drag 
coefficient  as  defined  by  Eqs.  (15)— (17)  of  Table  3.  Note 


that  Eqs.  (15)— (17)  of  Table  3  were  verified  with  water  bal¬ 
ance  measurement  [7].  Note  also  that  Eq.  (17)  approximates 
the  data  of  Springer  and  coworkers  [23]  since  X  =  1 1  @ 
ftd  =  0.9.  Also  note  that  Eq.  (17)  and  there  reference  [23] 
differ  from  the  data  presented  by  Fuller  and  Newman  [24] 
where  in  there  Fig.  4,  n&  (their  §)  is  constant  at  1.2  until 
X  =  4.  These  equations  will  vary  with  membrane  and  ionomer 
as  characterized  by  the  equivalent  weight  and  density  in 
Eq.  (19). 

The  species  transport  equations  (Eqs.  (3)-(6)  of  Table  2) 
are  solved  for  the  mass  flow  rates  of  the  hydrogen,  water,  and 
oxygen  species  based  on  the  mixture  component  velocities, 
u,  v,  and  w,  and  the  diffusion  mass  fluxes  J^\  for  every  time 
step  The  species  binary  diffusion  coefficients  are  calculated 
as  shown  by  Eq.  (14)  in  Table  3.  Fuller  and  Newman  [14] 
integrate  the  flux  expression  for  the  diffusion  of  water 
through  the  membrane  whereas  references  [13,15,17] 
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Table  1 

List  of  symbols 


Acv 

c 

CwK 

Dn,j 


Dw 

F 

I 

Io 


m n,K 

Mn 

mass„ 


nd 


w,K 

P 

Pn 

Q 

R 

S 


t 

fin 

T 

U,  V,  w 


Lcell 


X 

3f\y,K 


Activity  of  water  in  stream  K  (dimensionless) 

Specific  surface  area  of  the  control  volume  (cv)  (m- 1 ) 
Condensation  rate  (s-1) 

Concentration  of  water  vapor  at  K  interface  of  the  mem¬ 
brane  (molm-3) 

Binary  diffusion  coefficient  of  species  n  in  mixture  j 
(m2  s-1) 

Diffusion  coefficient  of  water  (m2  s-1) 

Faraday  constant  (96487C  mol-of-electrons-1) 

Local  current  density  (Am-2) 

Exchange  current  density  (Am-2) 

Mass  fraction  of  the  species  n  in  stream  K  (dimensionless) 
Equivalent  weight  of  a  dry  membrane  (kg  mol- 1 ) 
Molecular  weight  of  species  n  (kg  mol- 1 ) 

Mass  of  species  n  (kg) 

Electro-osmotic  drag  coefficient  (number  of  water 
molecules  carried  per  proton) 

Vapor  pressure  of  water  in  stream  K  (Pa) 

Pressure  (Pa) 

Partial  pressure  of  species  n  (pa) 

Volume  flow  rate  (m3/s) 

Universal  gas  constant  (8.314  J  mol-1  K-1) 

Source  term 
Time  (s) 

Membrane  thickness  (m) 

Temperature  (K) 

Velocities  in x,  y,  and  z  directions,  respectively  (ms-1) 
Cell  open-circuit  voltage  (V) 

Cell  voltage  (V) 

Channel  length  measured  from  anode  inlet  (m) 

Mole  fraction  of  water  in  stream  K 


Greek  symbols 
a 

h 

/ 

[i 

Pm,  dry 

P 

&  m 


Net  water  flux  per  proton  flux 
Permeability  in  the  £  direction 

Total  overpotential  for  oxygen  and  hydrogen  reactions 
(V) 

Water  content  in  the  membrane 
Dynamic  viscosity  (kg  s  m-2) 

Density  of  a  dry  membrane  (kg  m-3) 

Density  of  the  mixture  (kgm-3) 

Membrane  conductivity  (ohm- 1  m- 1 ) 


Subscripts  and  superscripts 


a 

c 

co2 

e 

H2 

K 

N2 

o2 

v 

w 

sat 


Anode 

Cathode 

Carbon  dioxide 

Electrochemical  reaction 

Hydrogen 

Anode  or  cathode  streams 

Nitrogen 

Oxygen 

Vapor 

Water 

Saturated 

Dummy  variable  for  direction  x,  y,  or  z 


assumed  a  linear  gradient  as  shown  by  Eq.  (15).  In  this  study, 
the  diffusion  coefficient  of  each  species  in  the  diffusion  layer 
was  reduced  arbitrarily  by  50%  to  account  for  the  effect  of 
porosity  and  pore-tortuosity.  The  flux  of  water  through  the 
membrane  is  critical  to  the  predictions  and  here  we  have 
used  the  same  equation  for  water  content  in  the  membrane 


as  given  in  Springer  et  al.  [11].  Eq.  (17)  gives  relationship 
between  the  electro-osmotic  drag  coefficient  and  water 
content  in  the  membrane.  This  equation  is  arrived  at  by 
curve-fitting  the  values  of  water  content  in  the  membrane 
and  electro-osmotic  drag  coefficient  presented  in  Springer  et 
al.  [1 1].  The  diffusion  coefficient  given  in  Eq.  (18)  in  Table  3 
is  taken  from  Nguyen  and  White  [15]. 

The  expressions  for  water  concentration  at  the  anode  and 
cathode  sides,  Cw?a  and  Cw,c>  are  calculated  according  to  Eq. 
(19),  and  the  activity  of  water  is  also  defined  in  Table  3.  More¬ 
over,  the  effects  of  local  hydrogen  and  oxygen  overpotentials 
are  accounted  as  shown  in  Eq.  (23).  It  is  important  to  note  that 
the  source  terms  in  Table  2  correspond  to  the  control  volume 
and  not  the  boundary  conditions  at  the  anode  or  cathode  inter¬ 
faces.  For  the  correct  determination  of  the  concentrations  and 
activities  at  the  membrane-diffusion  layer  interface,  the  mole 
fraction  for  each  species  used  in  these  equations  is  extrapo¬ 
lated  to  the  membrane  surface. 

It  is  important  to  note  that  here  the  transient  predictions 
below  assume  negligible  capacitance  per  unit  area.  That  is 
the  transient  current  is  the  sum  of  a  Faraday  current  density 
as  described  by  Eq.  (21)  of  Table  3  and  a  charging  current 
density,  7c(x,y).  This  charging  current  is  related  to  the  change 
in  overpotential,  rj (x,y): 


Ic(x,  y)  =  C 


y ) 

d  t 


2.1.  Numerical  procedure 

A  control  volume  technique  based  on  Star-CD  version  3.2 
commercial  flow  solver  was  used  to  solve  the  coupled  time 
dependent  governing  equations.  Since  we  neglect  7c(x,y)  the 
predictions  shown  below  would  apply  to  Pt  electrode  by  not 
a  Pt/Ru  electrode  since  Ru  has  a  large  C.  Several  subroutines 
were  added  to  calculate  and  account  for  the  source  terms, 
permeability,  electrochemical  reactions,  and  flux  of  protons 
and  water  across  the  membrane  according  to  the  equations  of 
Tables  2  and  3. 

Figs.  1  and  2  show  the  geometry  of  fuel  cell  modeled  in  this 
work,  which  consists  of  two  flow  channels  separated  by  diffu¬ 
sion  layers  and  MEA.  There  are  twenty  serpentine  passes  in 
the  flow  path,  so  that  the  flow  is  approximately  sixty  centime¬ 
ters  long  in  the  axial  direction  with  0. 1  cm(height)  x  0.08  cm 
(width)  cross-section  flow  area.  Each  diffusion  layer  has 
dimension  of  0.025  cm  (height)  x  3.20  cm  (width)  x  3.20  cm 
(length).  A  total  of  34  28  cells  x  200  28  cells  x  28  cells  (ele¬ 
ments)  were  used  to  model  the  fuel  cell.  That  is,  we  use 
34  cells  x  200  cells  x  6  cells  for  both  anode  and  cathode  dif¬ 
fusion  layers  and  34  cells  x  200  cells  x  6  cells  for  both  anode 
and  cathode  flow  channels,  34  cells  x  200  cells  x  2  cells  for 
the  MEA,  and  34  cells  x  200  cells  x  2  cells  for  caps  on  the 
flow  channels  corresponding  to  the  graphite  current  collec¬ 
tors. 

The  solution  procedure  used  in  this  commercial  flow 
solver  is  based  on  a  SIMPLE  algorithm  [25].  For  every  time 


Table  2 

Governing  equations  and  source  terms 


Governing  equations 


Conservation  of  mass 


Momentum  transport 


Mathematical  expressions 


dp  d(pu )  d(pv)  d(pw )  _ 

dt  +  dx  dy  dz 

d(pu)  d(pu)  d(pu) 

+  u -  +  v - 1-  w 


dt 

d(pv) 

dt 

d  (pw) 
dt 


dx 


dy 


d  (pv)  d(pv) 

+  u - +  v -  +  w 


dx 


dy 


d(pw)  d(pw) 

+  u - f-  v - b  w 


dx 


dy 


-~Sm 

d  (pu) 
dz 

d(pv)  _ 
dz 

d  (pw) 
dz 


dP  d  (  du  \  d  (  du\  d  (  du 

& +  &  Tto) +  ^  r^J +  ^  r& 1  +Spx 


dP  d  f  dv 
dy  +  dx  \  dx 


+ 


dv 
.  M — 

dy  V  dy 


d  (  dv 

+  to  Cto 


+  5 


py 


(1) 


(2) 


dP  d  (  dw\  d  (  dw\  d  (  dw\ 

to  +  to  (/*  to  )  +  to  C  to  )  +  &  V  to  J  +  Spz 


Non-zero  volumetric  source  terms  and  location  of  application  (see  Fig.  2) 

Sm  =  3h2  T  5aw  at  Z  =  Z3 
Sm  =  S O2  “b  *^cw  at  Z  —  Z2 


piv 


pLW 

77 


at  zi  <  z  <  Z4 


(7) 

(8) 


Hydrogen  transport 
(anode  side) 

Water  transport 
(anode  side) 

Oxygen  transport 
(cathode  side) 

Water  transport 
(cathode  side) 


d(pmu2 )  d(pmH  )  d (pmHo)  d(pmn2 )  d(JXt h2)  d(JVi h2)  d(JZtu2) 

+  u — —  +  v - 2 — —  +  w - - — —  =  — - — —  H - '±—  H - 7 - +  oh2 


dt 


dx 


dy 


dz 


dx 


dy 


dz 


d(pm.dw)  9(/7772aw)  9(/7772aw)  9(/7772aw)  d(Jx  aw)  9(*7y,aW)  d(Jz  aw) 

- I-  u - - - b  V - - - b  W - - -  =  - - - 1 - - - 1 - - - b  Saw 


dt 


dx 


dy 


dz 


dx 


dy 


dz 


d(pm  o2)  d(pmQ)  d(pm0  )  d(pm0 ,)  d(JXi  o2)  d(JVio2)  d(JZt  o2) 

+  U - - - b  v - 7 - b  w - - -  =  — - - 1 - - - f - - - b  Sq2 


dt 


dx 


dy 


dz 


dx 


dy 


dz 


3(/7772cw)  9(/7772cw)  3(/7772cw)  3(/7772cw)  d(Jx  cw)  d(JyCW)  9( -/y?cw ) 

- b  u - - - b  v - - - b  w - - -  =  — - - 1 - - - t - - - b  Scw 


dt 


dx 


dy 


dz 


dx 


dy 


dy 


(3) 


(4) 


(5) 


(6) 


I(x,  y) 

Sh2  = - mh2Acv  at  z  =  Z3 


a{x,  y) 

Saw  - - —  I(x,  y)MH2o  Acv  at  z  =  z3 


>o2 


i(x,  y) 
4  F 


Mq2  Acv  at  z  =  Z2 


1  +  2a(x,  y) 

S™  =  - — — —I(x,  y)MH2oAcv  at  z  =  z2 


2  F 


(9) 


(10) 


(11) 


(12) 
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Table  3 

Equations  for  modeling  electrochemical  effects 


Diffusion  mass  flux  of  species  /  in  £  direction 


-pDu 


dm kj 


(13) 


Binary  diffusion  coefficient  [26] 


2.334 


PDij(x,  y ) 


(• Pc-i  x  Pc-jPHTc-i  X  Tc-,fln((\IMi)  +  {\IMj))ll2 


=  3.64  x  10' 


Tcell 


y/ Tc-i  x  Tc-j 


(14) 


Net  water  transfer  coefficient  per  proton 


ofc  y)  -  n&(x,  y)-  tF  s  E>w(v  y)  ^Cwc^’  ^  Cwa^’  (15) 


I(x,  y ) 


m 


A.  =  0.043  +  17.8 l«a  —  39.85 al  +  36.0a;?;  0  <  aa  <  1 


Water  content  in  the  membrane 


(16) 


Electro-osmotic  drag  coefficient 


A  =  14  +  1.4(aa  —  1);  1  <  aa  <  3 


«d  =  0.0029AZ  +  0.05A  -  3.4  x  10“19  (17) 


Water  diffusion  coefficient 


Dw  =  nd5.5  x  10  11  exp 


2416 


1  1 


303  T 


(18) 


Water  concentration  for  anode  and  cathode 
surfaces  of  the  MEA 


Water  activity 


CvlK(x,y)  =  Pm’dry">^  '  n®~-  “  n-3  '• 


M, 


(0.043  +  17.8<3k  —  39.8aK  +  36.0<%);  a^  <  1 


m,dry 

Cwk(x,  a)  =  f^’dry  (14  +  1.4(<2k  —  1));  for  ak  >  1>  where  K  =  a  or  c 


(19) 


«K  = 


47m,dry 

^w,k^,  y) 

psat 

w.K 


(20) 


Local  current  density 


I(x,  y)  =  ffmfc  yVoc  -  Vceii  -  v(x,  y)l  (21) 


‘m 


Local  membrane  conductivity 


<ym(x ,  y)  =  (  0.00514^^CWa(^.  y)  -  0.00326  j  exp  (l268  ^  -  i  )  )  x  10z  (22) 


Local  overpotential 


RT 

n(x,  y)  —  -  In 

1  0.5  F 


I(x,  y)P(x,  y) 
l_7o0  2Po2(x,y) 


RT 

H - In 

0.5  F 


I(x,  y)P(x,  y) 
70h,  PH2 


(23) 


Viscosity  of  mixture 


/x 


=  m//x/  (24) 


step,  three  momentum  equations  corresponding  to  three  spa¬ 
tial  coordinates  are  solved,  followed  by  a  pressure  correction 
equation  that  corrects  the  mass  balance.  Species  transport 
equations  are  solved  after  the  bulk  flow  calculation.  The  mix¬ 
ture  properties  at  each  control  volume  are  calculated  based  on 
the  local  species  content.  The  anode  side  gas  mixture  contains 
hydrogen  and  water  vapor.  On  the  other  hand,  the  cathode 
side  gas  mixture  contains  oxygen,  water  vapor,  and  nitrogen. 
Therefore,  the  density  and  viscosity  of  the  two  flow  channels 
are  different  and  vary  from  one  location  to  the  other.  Note  that 
the  solution  procedure  for  the  time-dependent  predictions  is 
a  fully  implicit  scheme  [25].  Therefore,  the  newly  calculated 
values  of  the  variables  prevail  throughout  the  time  step  At. 
We  tested  the  results  to  be  sure  that  they  were  independent  of 
time  step  size.  That  is  for  example  we  stepped  from  0.02  to 
0.04  s  first  with  a  At  of  0.001  s  and  then  with  a  At  of  0.005  s 


and  the  results  were  the  same  to  within  5%.  Also,  during 
steady  state  simulations,  a  separate  grid  independence  test 
was  performed  by  increasing  and  decreasing  the  number  of 
the  grid  cells  on  a  straight  channel.  The  number  of  grid  cell 
was  decreased  and  increased  by  50%  of  the  base  case,  and 
predicted  results  were  compared  with  the  base  result.  The 
results  were  less  than  2%  different  from  each  other. 


3.  Results  and  discussion 

The  effects  of  the  rate  of  change  in  cell  voltage  on  cell 
performance  are  studied  for  three  different  rates  shown  in 
Fig.  3.  Condition  #1  is  a  step  change,  and  conditions  #2  and 
#3  correspond  to  incremental  changes  with  the  time  profiles 
detailed  in  Table  4.  For  most  of  predictions,  operating  flow 
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Time  (s) 

Fig.  3.  Cell  voltage  changes  with  time  used  in  this  study.  Condition  1  (•), 
a  step  change;  condition  2  (OX  an  increment  change  1;  condition  3  (V),  a 
slower  increment  change. 


Table  5 

Inlet  conditions,  physical  properties,  and  model  parameters 


Cell  voltage 

Volts  0.5 

Anode  channel  inlet  conditions  for  stoich  of  1.2 

Velocity  (ms-1)  2.6 

Mole  fraction  of  H2  0.729 

Mole  fraction  of  H2O  0.27 1 

Dew  point  temperature  (°C)  ~65 

Cathode  channel  inlet  conditions  for  stoich  of  2.0 

Velocity  (ms-1)  9.9 

Mole  fraction  of  O2  0. 165 

Mole  fraction  of  N2  0.622 

Mole  fraction  of  H2  O  0.213 

Dew  point  temperature  (°C)  ~57 

Operating  conditions 

Operating  pressure  (kPa)  101 

Cell  voltage  (V)  0.5 

Permeability  of  diffusion  layer  (m2)  3.3  x  10-15 


Time  (s) 

Fig.  4.  Transient  response  of  average  current  density  between  voltage  step 
change  and  voltage  increment  change  from  0.7  to  0.5  V. 

rates  for  this  study  corresponded  to  flow  rates  that  were  1.2 
times  greater  than  that  required  (by  the  measured  current)  for 
hydrogen  (i.e.,  20%  excess  hydrogen)  and  2.0  time  greater 
than  that  required  for  air  (i.e.,  100%  excess  air)  at  the  cell  volt¬ 
age  of  0.5  V.  Thus,  we  call  this  a  stoichiometric  flow  (stoic) 
of  1. 2/2.0  (i.e.,  1.2  for  anode  and  2.0  for  the  cathode).  The 
anode  and  cathode  flows  were  co-current  and  constant.  The 
initial  conditions  correspond  to  steady  state  concentrations 
and  velocity  and  pressure  distributions  at  0.5  V  with  a  feed 
velocity  of  2.6  m/s.  This  corresponds  to  an  initial  stoic  of 
2. 6/4.4  at  0.7  V  and  a  final  stoic  of  ~1. 2/2.0  at  0.5  V.  Thus, 


Membrane  thickness  (pan)  50 

Density  of  dry  membrane  (kg  m-3)  2000 

Equivalent  weight  of  a  dry  membrane  (kg  mol- 1 )  1.1 

Oxygen  exchange  current  density  (Am-2)  140 

Hydrogen  exchange  current  density  (Am-2)  1000 


Fig.  5.  Transient  response  of  average  current  density  with  different  voltage 
increment  change  from  0.7  to  0.5  V  ((•)  condition  #  2  current  density;  (O) 
condition  #  2  voltage;  (T)  condition  #  3  current  density;  (V)  condition  #  3 
voltage). 

we  have  an  excess  amount  of  fuel  and  air  at  time  equal  zero 
and  a  stoic  of  1. 2/2.0  at  time  equal  to  infinity.  Other  inlet 
conditions  and  physical  properties  are  shown  in  Table  5.  It 
is  important  to  note  that  these  conditions,  the  membrane  is 
close  to  saturation  (i.e.,  X  =  10)  during  all  of  the  transient 


Table  4 


Voltage  changes  profiles 


Condition  #1 

Time  (s) 

-0.10 

0.00 

0.01 

0.05 

0.10 

0.20 

0.30 

0.40 

0.50 

0.60 

0.70 

0.80 

0.90 

1.00 

1.10 

Voltage  (V) 

0.70 

0.70 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

Condition  #2 

Time  (s) 

-0.10 

0.00 

0.07 

0.08 

0.09 

0.10 

0.11 

0.12 

0.13 

0.14 

0.15 

0.16 

0.17 

0.18 

1.10 

Voltage  (V) 

0.70 

0.70 

0.70 

0.69 

0.68 

0.67 

0.66 

0.63 

0.60 

0.57 

0.54 

0.52 

0.50 

0.50 

0.50 

Condition  #3 

Time  (s) 

-0.1 

0 

0.05 

0.1 

0.15 

0.2 

0.25 

0.35 

0.45 

0.55 

0.7 

0.8 

0.9 

1 

1.1 

Voltage  (V) 

0.7 

0.7 

0.695 

0.69 

0.68 

0.67 

0.66 

0.63 

0.6 

0.57 

0.54 

0.52 

0.51 

0.5 

0.5 

Voltage  (V) 
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inlet 


outlet 


t  =  0.0s,  0.7V 


inlet 


outlet 


0.005  0.01  0.015  0.02 

Channel  width  (m) 

t  =  0.01s,  0.5V 


0.025  0.03 


inlet 


outlet 


t  =  0.08s,  0.5V 

Fig.  6.  Local  transient  current  density  (A/cm2)  contours  at  different  times  and  cell  voltages  for  a  step  voltage  change  condition  (i.e.,  condition  #1  of  Fig.  3). 


operation.  This  would  not  be  the  case  if  one  is  considering 
analysis  of  a  start-up  condition  with  a  dry  membrane.  Thus, 
the  fact  that  the  time  dependent  for  water  in  the  membrane 
is  approximately  zero.  This  as  discussed  below  that  transient 
response  is  not  a  function  of  the  changes  in  the  membrane 
hydration.  Note  that  our  experiment  for  PRIME  A®  ME  As 
indicates  that  the  dew  points  of  the  anode  and  cathode  are 
similar.  The  operating  pressure  is  101  kPa  and  cell  temper¬ 
ature  is  constant  at  70  °C.  The  membrane  thickness  used  in 
this  simulation  is  50  [xm. 

Fig.  4  shows  the  response  in  the  average  current  density 
when  the  changes  in  cell  voltages  are  step  and  incremental  as 


shown  in  Fig.  3.  The  cell  voltage  incremental  change  in  this 
figure  is  for  condition  #2  of  Fig.  3.  Results  are  shown  for  two 
anode  flow  rates  for  the  cell  voltage  change  of  condition  #2, 
corresponding  to  stoics  of  1.2  and  2.4  at  0.5  V.  In  this  figure 
when  the  cell  voltage  is  changed  in  one  time  step  from  0.7  to 
0.5  V,  the  average  current  density  rapidly  increases  from  0.35 
to  1.15  A/cm2.  Then  the  current  density  decreases  with  time 
and  reaches  steady  state  at  /  =  0.17  s  with  the  current  density 
of  0.90  A/cm2.  When  the  rate  of  voltage  change  is  slower 
as  for  condition  #2  with  an  anode  stoic  of  1.2,  the  current 
density  increases  initially  and  reaches  the  maximum  average 
current  density  of  1 . 10  A/cm2  at  /  =  0. 1 8  s.  Again,  the  current 
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inlet  outlet  inlet  outlet 


inlet  outlet  inlet  outlet 


Channel  width  (m)  Channel  width  (m) 

t  =  0.08s,  0.50V  t  =  1 0.0s,  0.50V 

Fig.  7.  Local  transient  oxygen  mole  fraction  contours  at  different  times  and  cell  voltage  for  step  voltage  change  condition  #1. 


density  decreases  until  t  =  0.35  s  when  it  reaches  the  steady 
state  value  of  0.90  A/cm2.  When  the  anode  stoic  is  increased 
to  2.4,  the  maximum  current  reaches  1.21  A/cm 2  at  t  =  0.18  s 
and  the  final  steady  state  is  0.92A/cm2.  Fig.  4  shows  that 
the  average  current  density  overshoots  its  final  steady  state 
value  for  conditions  #1  and  #2  and  that  the  peak  of  current 
overshoot  depends  on  the  rate  of  change  of  the  voltage  and 
the  amount  of  gas  flow  rate.  Higher  gas  flow  rates  result  in 
greater  current  overshoot  and  a  longer  time  to  reach  steady 
state. 

Fig.  5  illustrates  the  transient  current  response  with  the 
voltage  changes  of  conditions  #2  and  #3  of  Fig.  3.  This  figure 


confirms  that  the  magnitude  of  current  overshoot  is  affected 
by  the  rate  of  change  of  cell  voltage  and  it  is  consistent  with 
the  data  of  reference  [  1] .  Note  that  it  is  difficult  to  experimen¬ 
tal  by  obtaining  1  V/s  change  and  hence  the  model  provides 
useful  insight  to  extreme  conditions.  For  the  rate  of  change  of 
cell  voltage  corresponding  to  condition  #3  in  Fig.  3,  the  max¬ 
imum  average  current  density  reaches  0.92  A/cm2  at  t-  0.8  s 
and  cell  voltage  of  0.52  V.  The  current  density  reached  steady 
state  value  of  0.90  A/cm2  at  t  =  1 .0  s.  Additional  computations 
not  shown  here  verified  that  the  overshoot  could  be  eliminated 
completely  if  a  proper  function  for  the  cell  voltage  change  is 
implemented.  This  function  can  have  a  high  rate  of  voltage 
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Fig.  8.  Local  transient  current  density  (A/cm2)  contours  at  different  times  and  cell  voltages  for  increment  voltage  change  condition  #2.  Here  the  maximum 
current  occurs  at  /  =  0. 1 8  s. 


change  initially  as  long  as  the  voltage  change  is  small  as  the 
voltage  approaches  0.5  V.  The  reasons  for  the  overshoot  are 
discussed  below. 

Fig.  6  shows  the  local  current  density  on  the  membrane 
surface  for  four  time  steps  and  cell  voltages  for  the  step 
change.  These  figures  show  the  spatial  variations  of  local 
current  density,  which  cannot  be  provided  in  ID  and  2D  mod¬ 
els.  At  the  steady  state  corresponding  to  the  cell  voltage  of 
0.7  V  and  f  =  0.0s,  the  local  current  density  is  higher  in  the 
inlet  region  and  decreases  along  the  flow  channel  towards 
the  outlet  due  to  the  lowering  of  anode  activity  of  water  by 
electroosmotic  transport.  The  maximum  and  minimum  local 


densities  are  0.43  and  0.30  A/cm2,  respectively.  The  varia¬ 
tion  of  current  density  is  relatively  small  over  the  entire  cell 
area.  There  are  also  current  density  variations  at  the  channel 
bends  and  the  higher  current  density  is  larger  on  the  outside 
edges  of  the  bends.  This  is  because  the  sharp  turns  of  the 
bends  produce  re-circulation  zones  around  the  outward  cor¬ 
ners  (edges),  resulting  in  larger  velocities  and  hence  higher 
gas  concentrations  at  the  inner  corners  than  the  outer  cor¬ 
ners.  When  the  cell  voltage  is  lowered  to  0.5  V  in  0.01  s,  the 
local  current  density  significantly  increases  over  the  mem¬ 
brane  surface  but  the  contour  pattern  is  similar  to  the  /  =  0.0  s 
result.  The  highest  local  current  density  is  1.624  A/cm2  and 
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Fig.  9.  Local  transient  current  density  (A/cm2)  contours  at  different  times  and  cell  voltages  for  increment  voltage  change  condition  #3.  Here  the  maximum 
current  occurs  at  t  =  0.8  s. 


the  lowest  value  is  0.970  A/cm2.  As  time  increases,  there  is 
more  non-uniformity  of  local  current  density  distribution.  At 
0.08  s,  the  local  current  reduces  considerably  from  the  center 
of  the  membrane  surface  toward  the  outlet  due  to  decreasing 
of  oxygen  concentration.  The  highest  local  current  density 
is  located  around  the  inlet  with  a  value  of  1.74  A/cm2  and 
the  lowest  local  current  density  is  at  the  outlet  with  a  value 
of  0.63  A/cm2.  The  figure  for  10.0s  shows  that  at  the  final 
state  the  local  current  density  varies  from  1.71  to  0.27  A/cm2 
and  that  it  decreases  from  the  center  region  of  the  membrane 
surface  towards  the  outlet  region.  This  is  because  the  high 


reaction  rate  regions  consume  the  hydrogen  and  oxygen  at 
the  beginning  of  the  channel  and  the  down  stream  regions  are 
depleted  in  reactants.  Thus,  the  partial  pressure  of  O2  and  H2 
change  down  the  channel  and  under  ribs  even  though  X  and 
a  are  approximately  constant.  Again  for  other  operating  con¬ 
ditions  especially  those  from  a  drier  membrane  would  have 
different  profiles  and  the  model  may  need  an  accumulation 
term  for  the  water  balance.  Here  these  are  not  necessary  to 
illustrate  the  response  to  changes  in  the  load. 

Fig.  7  shows  the  transient  local  oxygen  mole  fraction  con¬ 
tours  on  the  membrane  surface  corresponding  to  the  current 


366 


S.  Shimpalee  et  al.  /  Journal  of  Power  Sources  156  (2006)  355-368 


inlet 


outlet 


inlet 


outlet 


inlet  outlet  inlet  outlet 


Fig.  10.  Local  transient  oxygen  mole  fraction  contours  at  different  times  and  cell  voltages  for  increment  voltage  change  condition  #2. 


density  of  Fig.  6.  These  figures  reveal  that  the  local  oxygen 
mole  fraction  patterns  are  consistent  with  the  local  current 
density  contours.  When  the  cell  voltage  is  changed  to  0.5  V 
with  a  time  step  of  0.01  s,  the  distribution  of  oxygen  mole 
fraction  becomes  quite  non-uniform  and  shows  a  rapid  deple¬ 
tion  towards  the  outlet  regions  of  the  membrane  surface.  The 
decrease  of  oxygen  mole  fraction  in  the  inlet  region  of  the 
membrane  cell  is  about  18%,  whereas  the  decrease  is  about 
30%  for  the  outlet  region.  Comparison  of  the  four  time  step 
contours  shows  that  the  oxygen  mole  fraction  in  the  inlet 
area  of  the  membrane  surface  reaches  its  steady  state  value 
of  0.09  (a  38%  reduction)  in  a  time  span  of  0.08  s  and  that 
it  reaches  the  steady  value  of  0.002  (a  98%  reduction)  over 


the  last  50%  of  the  membrane  in  10  s.  Under  these  condi¬ 
tions,  there  is  not  enough  oxygen  for  reactions  in  the  last 
section  of  the  cell  even  though  the  over  all  stoich  is  2.0  for 
the  cathode.  It  may  be  important  to  note  that  the  gradient 
of  O2  in  the  z  direction  (i.e.,  from  the  flow  channel  to  the 
electrode  surface)  are  much  smaller  than  the  changes  from 
inlet  to  outlet  for  these  value  of  the  diffusion  layer  shown  in 
Table  5. 

There  is  a  marked  difference  in  transient  local  current 
density  contours  for  the  incremental  change  of  cell  voltage 
corresponding  to  conditions  #2  and  #3  of  Fig.  3.  These  differ¬ 
ences  are  discussed  by  comparing  Figs.  8  and  9  with  Fig.  6. 
At  time  t-  0.0  s  the  current  density  contour  plot  is  the  same 
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Fig.  11.  Local  transient  oxygen  mole  fraction  contours  at  different  times  and  cell  voltages  for  increment  voltage  change  condition  #3. 


for  all  three  conditions.  Then  for  each  time  step  each  fig¬ 
ure  shows,  in  general,  that  the  current  density  contours  for 
gradually  decrease  from  the  inlet  to  the  outlet.  The  major 
difference  for  the  three  conditions  occurs  at  the  time  step 
at  which  the  maximum  overshoot  occurs.  The  magnitude  of 
the  current  density  for  all  three  conditions  is  comparable  in 
the  immediate  inlet  region  of  the  cell  but  a  slower  rate  of 
voltage  changes  allows  the  current  density  towards  the  out¬ 
let  region  to  approach  the  final  value  at  the  time  of  the  peak 
current.  At  this  peak  current,  the  current  density  distribution 
is  more  uniform  for  conditions  with  more  rapid  cell  voltage 
change.  The  current  density  distribution  is  most  non-uniform 


at  the  peak  current  for  condition  #3  (i.e.,  the  slower  rate  of 
cell  voltage  change).  As  a  result,  the  amount  of  overshoot  for 
the  average  current  density  is  lower  for  condition  #3.  These 
local  current  densities  are  a  result  of  the  concentration  of 
the  reacting  gases,  particularly  oxygen  under  these  operating 
conditions.  The  final  contours  of  both  Figs.  8  and  9  represent 
the  local  current  density  contours  at  the  steady  state.  These 
contours  are  similar  to  the  condition  #1  (i.e.,  a  step  change  in 
voltage). 

Figs.  10  and  11  show  the  transient  local  oxygen  mole 
fraction  for  incremental  voltage  change  of  conditions  #2 
and  #3,  respectively.  The  time  steps  and  cell  voltages 
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shown  in  Figs.  10  and  11  correspond  to  the  currents  of 
Figs.  8  and  9,  respectively.  These  figures  verify  that  the 
variation  of  local  current  density  for  the  specified  stoichiom¬ 
etry,  cell  voltage  and  time  step  change  is  due  to  the  local 
oxygen  concentration.  The  oxygen  decreases  from  the  cen¬ 
ter  of  membrane  surface  and  is  almost  depleted  toward  the 
outlet  when  cell  voltage  is  decreased  to  its  steady  state 
value. 

4.  Conclusions 

A  3D  time-dependent  simulation  of  a  PEMFC  was  shown 
for  different  rates  of  load  change.  Contours  of  local  current 
density  and  oxygen  concentration  were  presented  and  dis¬ 
cussed.  For  the  particular  operating  conditions  and  properties 
used  in  this  study,  the  membrane  is  close  to  complete  hydra¬ 
tion  and  the  results  indicate  that  there  can  be  current  overshoot 
of  the  final  cell  current  with  the  cell  voltage  is  changed  at 
1  V/s  from  0.7  to  0.5  V.  The  peak  of  this  overshoot  depends 
on  the  rate  of  voltage  change  and  anode  stoichiometry  and 
the  peak  of  current  overshoot  can  be  reduced  or  eliminated 
by  decreasing  the  rate  at  which  the  cell  voltage  is  changed. 
This  overshoot  behavior  is  a  result  of  local  non-uniformity  in 
the  current  density  distributions.  The  non-uniformity  of  local 
current  density  distribution  depends  on  anode  water  activity 
and  local  concentration  of  the  reacting  gases,  especially  the 
oxygen  concentration  for  the  conditions  of  this  study.  These 
results  can  help  our  on-going  studies  of  flow-field  configura¬ 
tion.  Design  changes  in  these  flow-fields  should  improve  the 
undesirable  fuel  cell  performance  such  as  the  current  over¬ 
shoot  as  a  result  of  load  change,  and  will  increase  the  perfor¬ 
mance  dramatically.  Also  simpler  models  that  approximate 
the  results  may  now  be  developed  and  compared  with  these 
results. 
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